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Abstract 

We compute the spectrum of normalizable fermion bound states in 
a Schwarzschild black hole background. The eigenstates have complex 
1/^ ■ energies. The real part of the energies, for small couplings, closely follow 

' a hydrogen-like spectrum. The imaginary parts give decay times for the 

, various states, due to the absorption properties of the hole, with states 

(N . closer to the hole having shorter half-lives. As the coupling increases, the 

}-j I spectrum departs from that of the hydrogen atom, as states close to the 

^ horizon become unfavourable. Beyond a certain coupling the lS'1/2 state 

is no longer the ground state, which shifts to the 2P3/2 state, and then to 
states of successively greater angular momentum. For each positive energy 
state a negative energy counterpart exists, with opposite sign of its real 
energy, and the same decay factor. It follows that the Dirac sea of negative 
^ , energy states is decaying, which may provide a physical contribution to 

• Hawking radiation. 

o\ ■ 

O ■ PACS numbers: 03.65. Ge, 04.70.Bw. 04.62. -|-v, 03.65.Pm 

O . 

^ ■ 1 Introduction 

O . Quantum theory in a black hole background has been extensively studied by 

many authors. Detailed discussions of this problem are contained in the books 
5-H . by Birrell & Davies and Chandrasekhar ,2, , and the review paper by Brout et 

^f^' al. in]. Much of the attention in this work is focussed on the wave equation and 

its scattering properties. Detailed studies of the Dirac equation in a black hole 
background are less common. Indeed, the lowest order scattering cross section 
for a fermion in a black hole background has only recently been computed 
' In this paper we investigate another previously neglected aspect of quantum 

mechanics in a black hole background. This is the existence of the bound state 
spectrum for particles orbiting a spherically-symmetric point source. That is, 
we study the gravitational analogue of the hydrogen atom orbitals. 

There has been strangely little effort devoted to the study of the bound state 
spectrum, despite the fundamental importance of the electromagnetic analogue. 
But it is clear that these states must exist — how else can one provide a quantum 
description of a particle in orbit around a black hole? These states must also 
be essential in the quantum description of the capture process. The problem 
was discussed in 1974 by Deruelle and RufRni jBj, who described the existence 
of resonance states in the Klein-Gordon equation. Further significant contribu- 
tions were made in a series of papers by Gaina and coauthors [3IHlin]. These 
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papers give various analytic expressions for the real and imaginary parts of the 
energy in a series of limiting cases. 

Much of the study of quantum mechanics in a black hole background has fo- 
cussed on the related, though distinct, problem of finding the quasi-normal mode 
spectrum. Quasi-normal modes are purely ingoing at the horizon, and outgoing 
at infinity |2j. These boundary conditions produce a spectrum of eigenstates 
with complex-valued energies. The significance of these quasi-normal modes 
comes from their use in describing black hole oscillations. But the boundary 
condition at infinity implies that these modes are not normalizable, so cannot 
represent bound states. The problem of interest here is to find these bound 
states, so we seek solutions which are purely ingoing at the horizon, and which 
fall off exponentially at infinity. 

For a particle of mass m in the field of a black hole of mass M the dimen- 
sionless coupling strength is defined by 

mM 

a= — - 1) 

where rup is the Planck mass. In this paper we compute the fermion bound 
state spectra for a in the range • • ■ 6. If the bound particle is assumed to be an 
electron, this range corresponds to black holes of masses up to 1 x lO^^kg, which 
is the scale appropriate for primordial black holes formed in the early universe. 
Computing the energy spectrum is more complicated than the hydrogen atom 
case for two main reasons. The first is that the radially-separated Dirac equation 
contains three singular points, only two of which are regular. There is no special 
function theory appropriate for the study of such equations, so we have to resort 
to a range of numerical techniques to find the spectrum. The second problem is 
that the singularity at the centre of a Schwarzschild black hole acts as a current 
sink. All normalizable states must therefore decay in time, and we must search 
for eigenstates over the two-dimensional space of complex energies. The states 
we construct therefore all have a finite half-life, so can be viewed as resonance 
states. The interpretation of these states is discussed in sectional 

Despite these difficulties, the problem can be tackled numerically, and we 
present a range of results for the real and imaginary parts of the energy. These 
are sufficient to predict how the spectrum will behave for larger values of the 
coupling constant. The first result, which is entirely to be expected, is that the 
orbitals become increasingly tightly bound as the coupling increases. It follows 
that, for a given state, the energy will initially decrease with a, but will eventu- 
ally turn round and start increasing as the particle spends too much time inside 
the classical radius of minimum energy. States with higher angular momentum 
then become energetically favourable as a increases. For example, we show that 
beyond a w 0.6 the 15*1/2 state is no longer the ground state. While the real 
part of the energy behaves in quite a complicated fashion, the imaginary part, 
which controls the decay rate, simply increases in magnitude. This is also as one 
would expect. The closer the orbital density is to the singularity, the greater 
the probability of capture. 

We start by discussing the Dirac equation in a Schwarzschild background in 
an arbitrary gauge. This is helpful in establishing a range of gauge-invariant 
results. In particular, the energy conjugate to time translation symmetry is con- 
firmed to be a gauge invariant quantity. This is important in order to guarantee 
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that the quantity is a physical observable. We next establish the behaviour of 

the wavefunction around the horizon, which is sufficient to establish that the 
states must decay exponentially with time. We then turn to a specific choice 
of gauge that is well-suited to numerical solution. We solve the equations by 
simultaneously integrating out from the horizon and in from infinity. We then 
vary the energy to ensure that the solutions match at some finite radius. This 
process guarantees that we find a global, normalizable bound state. A set of 
spectra are obtained, and the density is plotted for a range of states. Decay 
rates and expectation values for the distance are also presented. We end with 
a discussion of the significance of these bound states, and the possible physical 
processes that they may generate. Except where stated otherwise, natural units 
with G = h = c = 1 are assumed throughout. We employ a spacetime metric 
with signature —2. 



2 The Dirac equation 

We start by defining a general parameterisation of the Schwarzschild solution. 
This general form will help to guarantee that various expressions are gauge 
invariant. We let {70571,72573} denote the standard gamma matrices in the 
Dirac-Pauli representation, and introduce polar coordinates {r,6,(j)). Prom 
these we define the unit polar matrices 

7r = sin^(cos07i + sin(^72) + cos^73 
70 = cos0(cos07i + sin(/)72) — sin^73 
70 = - sin^ 71 + cos^ 72 . (2) 



In terms of these we define the four matrices 



5* = oi7o - a2lr 9 =- -le 



g'' = -&i7r + ^270 5'^ = T-570 (3) 

where (ai, 02, &i, ^2) are scalar functions of r satisfying 

aib\ — 0262 = 1 

{Wf - (62)' = 1 - 2M/r. (4) 
The reciprocal set of matrices are therefore 

gt = ^i7o - ^^27^ ge = r^0 

gr = aijr - a27o 94, = rsme^4,. (5) 



These matrices satisfy 



{5^5^} = 2.g''''/ 

{9ti,gi^} = '^9(^1^1 

{5^54 = 25^;/ (6) 
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where /i, run over the set (t,r,9,(f>), I is the identity matrix, and g^i, is the 
spacetinic metric. The hne element defined by this metric is 

g^t.dx^'dx'' = (l - 2M/r)df + 2(0162 - a2bi)dtdr - {{aif - {a2f)dr^ 

-r^{d0^ + siTi^0d(l)^). (7) 

This hne element is the most general form one can adopt for the Schwarzschild 
solution. There is only one degree of freedom in equation jT)), since the terms 
are related by 

(1 - 2A//r)((ai)2 - (03)') + (0162 - a2hf = 1. (8) 

This arbitrary degree of freedom corresponds to the fact that the time coordinate 
is only defined up to an arbitrary radially-dependent term. That is, we can set 

i=t + a{r), (9) 

and the new line element will be independent of the new time coordinate t. 
Rather than think in terms of changing the time coordinate, however, it is 
simpler for our purposes to always label the time coordinate as t and instead 
redefine ai and 02. These then transform as 

01 I— > ai = ai ~ 620^' 

02 1-^ 02 = 02 - ha' , (10) 

with 61 and 62 unchanged. Throughout dashes denote derivatives with respect 
to r. It is straightforward to confirm that the new set (ai, 02, 61, 62) still satisfy 
the constraints of equation Q . 

The four variables ai, 02, 61 and 62 are subject to two constraint equations, 
so must contain two arbitrary degrees of freedom. The first arises from the 
freedom in the time coordinate as described in equation 1)10(1 . The second lies 
in the freedom to perform a radially-dependent boost, which transforms the 
variables according to 



0-1 bi\ /cosh/3 sinh/3\ / ai bi 
0-2 ^2/ lsinh/3 cosh/3y \a2 62 



(11) 



where (3 is an arbitrary, non-singular function of r. This boost does not alter 
the line element of equation Q. Outside the horizon we have |&i| > \b2\, and 
in the asymptotically fiat region bi can be brought to -|-1 by a suitable boost. 
It follows that we must have 

bi>0 Vr > 2M. (12) 

At the horizon we therefore have bi positive, and 62 — ±61. For black holes 
(as opposed to white holes) the negative sign is the correct one, as this choice 
guarantees that all particles fall in across the horizon in a finite proper time. 
This sign is also uniquely picked out by models in which the black hole is formed 
by a collapse process. We can therefore write 

62 = -61 at r = 2M. (13) 
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Combining this with the identity oifei — 0262 = 1 we find that, at the horizon, 
we must have 

0162 - 02^1 = -1 at r = 2M. (14) 

The diagonal form of the Schwarzschild metric sets 0162 — a26i — 0, so does 
not satisfy this criteria. But for this case the time coordinate t is only defined 
outside the horizon, and the horizon itself is not dealt with correctly. 

We now have a general parameterisation of the Schwarzschild solution in 
an arbitrary gauge. The next step is to write down the Dirac equation is this 
background. This is 

= wV', (15) 

where 

V^V = (9^ + ^rf S„0)V', I],^ = ^[7„7^]. (16) 

The components of the spin connection are found in the standard way (see 
Nakahara ^Oli for example). These turn out to give 



9^\vf,,^, = ( ^2 + ^) 70 - f + ) >. (17) 



For the Dirac spinor we use a radial separation of the form 



where E is the (complex) energy and 

CTr = sin0(cos(/) (Ti + sin0 (T2) + cos6' 0-3. (19) 

The angular eigenmodes are labeled by k, which is a positive or negative nonzero 
integer, and /i, which is the total angular momentum in the 9 = direction. 
Our convention for these eigenmodes is that 

i(T-L + h)x^^nnx^,., ^ = ...,-2,-1,1,2,.... (20) 

The positive and negative k modes are related by 

^.X^ = xIk- (21) 
The trial function (|18() results in the pair of coupled first-order equations 
h b,\ fu[\ ^ /uA ^22) 



(23) 



b2 bij V"2/ \"2 

where 

K/r ~b[/2 + ia2E i{m + aiE) - h'2/2' 

~i{m~ aiE) -b'2/2 -n/r - b[/2 + ia2E 

These are the equations we wish to solve for complex energy E. It is first worth- 
while confirming that the equations are gauge invariant. A redefinition of the 
time coordinate is equivalent to the transformations described in equation (|10|l . 
These are combined with the transformation 

Ml uie-'^°' U2 i-> U2e-'^" (24) 
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which together ensure that equation H22|) is still satisfied. The radial boost 
defined by equation is combined with the transformation 

/cosh(/3/2) -sinh(/3/2)\ /uA 
\~smh{f3/2) cosh(/3/2) J \u2j ^ ^ 

to again ensure that the equation is still satisfied. In either case we see that the 
eigenvalue E is unchanged, so is a true gauge-invariant quantity. 

The angular separation of equation H18() is clearly justified from the form of 
the Dirac equation. The separation into energy eigenstates is gauge invariant, 
but it is helpful to see the separation in a gauge where the Dirac equation takes 
on a Hamiltonian form. This is provided by the 'Newtonian' gauge which 
sets 

ai = 1 a2 — 

bi = l b2 = -(2M/r)i/2. (26) 

In this gauge the Dirac equation takes on the simple form 

^^^-^7'^[^y\l + l)^ = mi., (27) 

where ^ is the Dirac operator in flat Minkowski spacetime. This equation is 
manifestly separable in time, so has solutions which go as exp{—iEt). Since the 
separation works in this gauge, it must work in all others. We will return to 
this gauge choice when we turn to finding the energy spectrum. 

The nature of equation H22() can be understood more clearly by writing it in 
the form 

('-^"/••)(::0-(4 ^)«(::)' <^«> 

This exposes the fact that the horizon is a regular singular point of the radial 
equations. The same is true of the origin, and infinity turns out to be an irregular 
singular point. This implies that the radial equations cannot be manipulated 
into second-order hypergeometric form, as one is able to do for the hydrogen 
atom. The closest the equations come to a recognisable form is that of Heun's 
equation, which generalises the hypergeometric equation to the case of four 
regular singular points on the complex plane [r2) . But Heun's equation can 
usually only be analysed using numerical techniques, and these are the tools we 
will apply to equation (|22|l 

The presence of singular points means that we must check carefully that our 
solutions behave appropriately at these points. The point at infinity is not an 
issue, as we seek solutions that decay exponentially. Similarly, the origin is not 
a problem. We expect that the function will be weakly singular there as the 
origin acts as a current sink, and this is indeed the case. The horizon, however, 
is more complicated. The wavefunction must be well behaved at the horizon 
if it is to represent a physical solution. To test this we introduce the series 
expansion 

oo oo 

Ml = 77"^afc?7^ U2 = ri'^f3kv'', (29) 

fc=0 k=0 
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where rj = r — 2M. On substituting this series into equation (|28|) . and setting 
rj = 0, we obtain the indicial equation 



det 



= 0, (30) 

r=2M 



r 

where / is the identity matrix. Employing the result that 

bib[ - 6262 = ^I/r^ (31) 
we find that the two solutions of the indicial equation are 

s = 0, -i + AiME{bia2 - b2ai)r=2M- (32) 
Equation H14|) then tells us that the two indices are 

s = 0,-^ + AiME. (33) 

These indices are therefore gauge invariant. The regular root s = ensures that 
we can always construct a solution that is finite and continuous at the horizon. 
The singular branch gives rise to discontinuous solutions with an outgoing cur- 
rent at the horizon. These can be used to provide a heuristic explanation of the 
Hawking radiation jJT] . It is clear that the non-zero indicial root gives rise to a 
wavefunction that is ill-defined at the horizon, and so cannot represent a phys- 
ical state. We must therefore confine our search for bound states to solutions 
that are regular at the horizon. 

The regular and singular solutions are related by a generalized form of time- 
reversal symmetry. For this we define 

^{t, x) - ^—^^1^-^(61^0 - hjrWi-t + f{r), x), (34) 

which effectively reverses the time direction using the normalized timelike Killing 
vector. In terms of the ui and U2 functions, the new solution is characterised 

by 



Ul 



eM-^E*f{r)) f bi 62 



(35) 



^U2) (l-2A//r)i/2 \-b2 -bi^ 
where /(r) is determined by 

(1 - 2M/r)drf{r) = 2(ai&2 - a2&i). (36) 

The time-reversed solution has energy E* and so is exponentially growing in 
time. The solution is also is singular at the horizon, employing the non-zero 
root of the indicial equation, and is not normalizable. 

Eigenmodes with different values of k and fi are orthogonal. For states with 
the same values of k and fi the quantum inner product can be taken as 

/•OO 

= / dr {ai{ulvi + ulv2) + a2{ulvi + ulv2)) , (37) 
Jo 

where the Ui and Vi denote the radial functions in ip and 4> respectively. Current 
conservation for the Dirac equation is summarised in the relation 



= {biiuiu*2 + U2UI) + b2{uiul + M2"2)e"'^^"^*^*] • (38) 



{uml + M2U2) + a2{uiu*2 + U2ul)e-'''^ 
d 
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Again it is straightforward to confirm that this equation is gauge invariant. The 
right-hand side of this equation defines times the radial flux. We denote this 
by J, 

J(r) — hx{uxu*2 + U2U1) + &2(uiwJ + uiu^)- (39) 

For spatially normalizable states we must have J 1-^ as r i— > 00. But at the 
horizon we also have 

J = -6i|ui-u2p, (40) 

which defines an inward-pointing current. At the horizon, the regular solution 
has 

|«i-U2p = |ao-/3o|', (41) 
using the power series expansion of equation H29(l . The coefficients are related 

by 




It is therefore impossible to satisfy ao — /3o for finite energy, so there must be 
a non-vanishing inward current present at the horizon. This in turn tells us 
that the state must decay. This decay takes place at the origin, where unitary 
evolution breaks down 113) . For bound states the energy E must contain 
real and imaginary terms, so we set 

E = uj~w. (43) 

Current conservation now takes the form 

^ = 2i'(ai(uiul + U2M2) + 02(^1^2 + U2ul)). (44) 
ar 

Given a set (wi, U2, E, k) which solve the radial equation H22|l a new solution 
set is generated by the transformation 

(ui, U2, E, k) 1-^ {u2,uI,-E*,~k). (45) 

It follows that the real part of the energy spectrum is symmetric about the 
zero point. That is, for a state with real energy w a corresponding antiparticle 
state exists with real energy ~lu. The decay rate is the same for both states, 
however. If we assume that the vacuum is constructed from the Dirac sea of 
negative energy states, then this vacuum will decay in time. A loss of negative 
energy states can be equally interpreted as generation of positive energy states, 
which provides a suggestive physical model for Hawking radiation. 

3 The energy spectrum 

To solve for the energy spectrum we work mainly in the Newtonian gauge of 
equation (|26)l . In this gauge the interaction with the black hole is defined solely 
by an interaction Hamiltonian Hi given by 

(46) 
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Dimensional constants arc included in a number of equations in this section to 
illustrate certain features of the problem. The line element for the Newtonian 
gauge has flat spatial sections for constant t, so the quantum inner product 
between states has the simple flat-space form 

(^10) = J d^x^p^(j). (47) 

The interaction Hamiltonian is not Hermitian, as we have 

Hi - hI = -ih(2GMr^Yl^5{x). (48) 

It is straightforward to check that all wavefunctions approach the origin as 
r"'^/'*, so the non-Hermitian part of Hi has finite expectation. This confirms 
that Hermiticity only breaks down at the origin, as stated earlier. In this respect 
it may be more natural to refer to Hi as a 'pseudo-Hamiltonian', in the sense of 
an operator acting on an open quantum system ^1] . There is no doubt that the 
system described here is open, as the singularity is not treated as part of the 
quantum system. But the system is only open in an extremely simple fashion. 
There is no ambiguity in either the time evolution of the state, or the correct 
definition of the observable energy. Time evolution is defined by the Dirac 
equation, in whichever gauge one chooses to adopt, and the energy is defined by 
the energy- momentum tensor. The fact that we have a Hamiltonian description 
at all is a result of a series of gauge choices, so one must be careful not to place 
too strong an interpretation on the gauge-dependent quantity Hi. 

The bound state energy eigenspectrum is determined entirely by the prop- 
erties of the wavefunction at the horizon and at infinity. The demands that 
the wavefunction is finite at the horizon and falls off exponentially at infinity 
are sufficient to produce the spectrum. But it is only by considering the global 
properties of the wavefunction that the imaginary contribution to the energy is 
fully understood. Decay only takes place at the singularity, and the decay rate 
for a given eigenstate is naturally related to the behaviour of the wavefunction 
near the singularity. If the spatial degrees of freedom in an energy eigenstate 
are normalized such that 



dr {uiu\ -f U2ul) = 1 (49) 

then the imaginary component of the energy, ^iv, is determined by 

n(2GM)i/2 1 , 

^"r™ 2 ;:57j("i"i + "2W2)- (50) 

This identity only holds if the state is globally normalized. It provides a fur- 
ther independent check that the solutions we obtain numerically are globally 
normalizable bound states. The decay rate increases for states with a larger 
probability density near the singularity, as one would expect. The fact that the 
states approach the origin as r~^/^ ensures that the radial probability density 
tends smoothly to zero at the singularity. The presence of the singularity does 
not prevent the formation of normalizable states, and the singular nature of the 
wavefunction is no worse than that of the ground state of the hydrogen atom. 

The interaction Hamiltonian is independent of the speed of light, so the 
non-relativistic approximation to the Dirac equation results in the Schrodinger 
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equation 



V'NR + iri\ ^ I ^371^ V V^nrJ = -E'NrV'nr, (51) 



2m 



where the subscript NR denotes non-relativistic. If we now introduce the phase- 
transformed variable 



4- 



where 



we see that ^ satisfies 



V'NRexp(-i(8r/ao)i/2) (52) 



^2 



* * = £;nr*- (54) 



2m 



In the non-relativistic limit the energy spectrum is therefore given by the grav- 
itational analogue of the hydrogen atom spectrum [S], 

C^A'Pm^ 1 

-Enr = —3 J, 77 = 1,2,.... (55) 

In terms of the Planck mass rup we can also write 

The fact that we have a reasonable starting point for the spectrum in the weak- 
coupling limit is valuable, as our method involves searching for eigenvalues over 
the complex energy plane. By analogy with the hydrogen atom case, we expect 
that the non-relativistic spectrum will be a reasonable approximation provided 

Mm 

< 1. 57 

m^ 

Returning to the full, relativistic equation H22|l . we convert this to dimen- 
sionless form by introducing the dimensionless distance variable 



(58) 



which ensures that the horizon lies at a; = 2. We also introduce the dimensionless 
coupling coefficient 

a = — ^ (59) 
ml 



and energy 

EM 



c^mp 



(60) 



10 



In terms of these our eigenvalue problem becomes 



(x-2) 



1 (2/a;)i/2 
(2/.x)i/2 1 



K ix{a + e) ~ (8a;) ' \ / ui 

ix{a — e) — (Sx)^-^/^ — k J \U2 



(61) 



where the dashes now denote derivatives with respect to x. We seek eigenvalues 
e for fixed coupling a. 

Two complementary methods are employed to solve the eigenvalue problem. 
We start with a series expansion around the horizon of the regular branch of 
the solution. The restriction to this branch removes two degrees of freedom at 
the horizon, so the function is uniquely specified up to an overall magnitude 
and phase. These are chosen conveniently by setting ui = 1 at the horizon. 
The power series expansion extends the solution a short distance away from 
the horizon, from where the values of (mi, U2) can be used to initiate numerical 
integration of the differential equation (|61|l . For most values of e the numerical 
integrator will start to increase exponentially after a finite distance. The aim 
initially is to vary e so as to push this distance out as far as possible. This 
requires a reasonable initial guess for the eigenvalues, which is where the non- 
relativistic approximation is helpful to get things started. 

Once we have achieved a reasonably accurate value for e, we turn to a more 
sophisticated method to improve accuracy. We seek normalizable states for 
which Tp is finite over all space. To be confident we have found such a state we 
need to numerically integrate inwards from infinity, as well as outwards from the 
horizon. If the solutions for ui and U2 can be arranged to match at some suitable 
radius then we have found a global solution to the first-order equations (j61|l . To 
expand about infinity we need to take care of the essential singularity present 
there. A suitable series expansion is provided by 

where 

2 2 2 I 2 i j:i2\ /co^ 

p = a — e = — 1— r(?Ti c — )■ (63) 

The definition of p involves a complex square root, and the branch is chosen so 
that p has a positive real value, ensuring the wavefunction falls off exponentially. 

The fact that only one root of the indicial equation is used implies that, for 
a given tjj is specified at infinity up to an arbitrary magnitude and phase. The 
first few terms in the series expansion (|62(l are used to compute -0 at a finite 
radius and these values are then numerically integrated inwards. A certain 
amount of fine tuning is then required to pick the radius at which to attempt 
matching. Once a radius is chosen the matching condition is that the inward and 
outward values of the two complex functions ui and U2 agree. This condition 
is converted into a set of four scalar equations which state that the real and 
imaginary differences vanish. In addition we have four arbitrary parameters to 
vary — the real and imaginary terms in the energy, and the magnitude and 
phase of the function integrated inwards from infinity. This system of four 
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equations and four unknowns is then solved by a Newton-Raphson method. 
This converges very quickly and affords good control over accuracy. 

Three independent checks were performed on the energy spectrum achieved 
by this method. The first was that the calculations were repeated using the 
same scheme in a different gauge. The gauge chosen for comparison is defined 
by advanced Eddington-Finkelstein coordinates, with 

ai = l + M/r a2 = M/r 

61 = 1 - M/r 62 = -M/r. (64) 

The second test involved using a minimax routine to find the energy spectrum. 
This method is less accurate, but gave good agreement for the states of lowest 
energy. The final check was to confirm that, after normalization, the states 
satisfy the identity of equation H5()|l . This check was again satisfied to high 
precision. 



4 Results 

The real parts of the energy for the three lowest-energy states are plotted in 
figure n The vertical axis plots the real part of the energy in units of the rest 
energy of the particle, which is given by 

(65) 

The fact that we obtain this dimensionless ratio reflects the equivalence princi- 
ple. The mass m does not effect the spectrum on its own — the spectrum only 
depends on the product mM . States are labelled using the standard spectro- 
scopic scheme. In this scheme k = 1 corresponds to 5*1/2, = 2 to ^3/2 s-nd 
K = — 1 to Pi/2- For each eigenvalue k a ladder of levels is obtained. 

The energy spectrum illustrates a number of remarkable features. For small 
a the spectrum resembles that of a hydrogen atom. But as the coupling increases 
the energy of the lS'1/2 state reaches a minimum and then starts to increase. The 
gravitational case avoids the Z = 137 catastrophe of the relativistic Coulomb 
problem. This is to be expected — coupling strengths with a > 1 are routinely 
achieved astrophysically and such objects appear to be stable. We also see that 
as a increases beyond 0.6 the P3/2 state appears to take over as the ground 
state. This is confirmed in figure 13 which shows the spectra of the 5*1/2, ^3/2 
and states out to a = 1.4. We see clearly that around a = 0.6 the 2P3/2 
state takes over from 15*1/2 ^ the ground state, only to be replaced in turn by 
the 31)5/2 state at a = 1.2. An explanation of this phenomena can be found in 
the classical expression for the binding energy in a Schwarzschild potential. 

For a particle of mass m in a Schwarzschild background the dynamics reduces 
to motion in the effective radial potential 

GMm J2 / 2GM\ 



where J is the angular momentum of the particle. This is illustrated in figureOl 
For J > \/T2GMm/c, classical bound states can exist as the effective poten- 
tial has a minimum, but if the particle's angular momentum is smaller than 
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Figure 1: The real part of the bound state energy, in units of mc^. The horizon- 
tal axis labels the dimensionless coupling coefficient a, and the lines represent 
the value of the energy for the coupling at the left of the line, with a ranging 
from 0.1 to 0.6 in steps of 0.05. The S1/2, P1/2 and P3/2 orbits are shown. 
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Figure 2: The energy spectra of the 5'i/2, -P3/2 and D^j^ states. At around 
a = 0.6, the 2P state becomes the ground state, and beyond a = 1.2 it is 
replaced by the 3D state. 



'12GMm/c it becomes insufficient to support a classical orbit. For a circular 
orbit at radius r the conserved relativistic energy, conjugate to time translation, 
is 

2 r-2GM/c'^ 
^ = ^ %i/2(,_3GM/c^)i/2 - (67) 

The radius r and angular momentum J are related by 

J2_ _ GMr^ 

r- iGM/c" ■ ^ ^ 

Now suppose we attempt a form of naive Bohr quantisation by setting 

J = nU. (69) 
Converting to dimensionless quantities the effective potential becomes 

po) 

and the orbital energy is 

i = (71) 

" (a;(:E-3))'^' 



where 
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Figure 3: Bohr quantisation of classical circular orbits. The left hand plot shows 
the effective potential in units of mc^ for a — 0.5 and J ~ nh, with n — 1,2, 3. 
For n = 1 no classical bound state is possible. As n is increased a minimum 
forms in the potential, and the barrier between the minimum and the singularity 
grows larger. The right-hand plot shows the relativistic energy in units of mc^ 
as a function of a. The first eight n values are shown. For each n value the 
minimum energy is achieved when a = n/^/T2. 

In the small a regime this reproduces the spectrum of equation H5fi|l. But as a 
increases the energy falls to a minimum at = beyond which the orbit 

no longer exists for a given n (see figure EJ- The minimum energy achieved is 
0.94mc^, corresponding to a: = 6. Inside this radius no stable classical circular 
orbits exist. In the quantum description we find that as a increases the orbits 
get more tightly bound around the horizon. As the coupling increases the orbits 
are dominated by terms inside a: = 6 and so become energetically less favourable. 
The ground state is then one of higher angular momentum, for which the orbit 
is less tightly bound. Figure [21 also shows that as a increases the 151/2 state 
becomes unbound. This effect is also seen classically, as circular orbits with 
r < 4M are known to be unbound, as well as unstable. 

The form of the effective potential illustrates a further feature of the quantum 
states, which is that the quantum decay can be interpreted as a tunnelling 
phenomena. This is certainly a valid picture for states with n > ^/T2a. For a 
fixed n, as a increases, the potential barrier decreases and we expect that the 
tunneling rate onto the singularity will increase. This is indeed the case, as we 
discuss further in sectional 

Figure ^ shows how states of successively higher angular momentum take 
over as the ground state as the coupling is increased. In the small a/n limit, 
the energy levels resemble those of the classical orbits. At larger couplings, 
the energy of a given state falls to a minimum, and then begins to increase 
again, apparently without limit. The a value at which the minimum occurs is 
roughly proportional to the angular momentum of the state, with a = 0.58k — 
0.10 providing a good fit. The maximum binding energy available increases 
with angular momentum, to beyond E — 0.88mc^. This means that quantum 
mechanics predicts around twice the classical value for the radiation efficiency 
of accretion processes. To confirm this effect we need to find the limiting value 
of the binding energy for astrophysical values of a. For an electron around a 
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Figure 4: Energy levels of states with higher angular momenta. This plot shows 
the energy levels of the lowest energy states with a range of angular momenta 
K = 1 • • • 10. It illustrates how each state takes a turn as the groundstate, as the 
coupling is increased. The positions of the energy minima are linearly-spaced 
in a and have increasingly large binding energies. 

solar-mass black hole, for example, we have a = 4 x 10^^, so a large a limit of 
our equations should be very accurate. 

5 Wavefunction properties 

With our current choice of gauge the radial form of the wavefunction is best 
visualised by plotting times the timelike component of the current. We denote 
this p, so 

P=\ui\^ + \u2\^. (73) 

The gauge invariant definition of p is that it is times the density as measured 
by observers in radial free-fall from rest at infinity. The first four S1/2 states 
for small coupling are shown in figure The plots are very similar to those for 
the non-relativistic hydrogen atom. In all cases the peak of the wavefunction 
is a long way outside the horizon, with only a small fraction of the probability 
density lying inside the horizon. 

As we increase the coupling to a = 0.35 we obtain the series of plots in 
figureini Predictably, the wavefunctions start to bunch in closer to the horizon. 
Slightly more surprisingly, the nodal structure disappears for larger couplings. 
The density no longer drops down to near zero at a number of nodes, but 
instead a number of dips are present. If we increase the coupling further still, to 
a — 0.5, the dips themselves are largely washed out and we obtain the somewhat 
structure-less plots shown in figured 

Some additional insight into the nature of the orbitals is obtained by calcu- 
lating the expectation value of r. With our current gauge choices this is defined 



16 



1S 2S 




500 1000 500 1000 1500 2000 



3S 4S 




1000 2000 3000 4000 2000 4000 6000 



Figure 5: The radial probability density for the lSi/2, 25*1/2) 351/2 and 45i/2 
states for a coupling of a = 0.1 . The horizontal axis is the dimensionless radius 
X, and the horizon lies at a; = 2. All plots are started from the horizon. The 
part of the density inside the horizon is not plotted, though in all cases this 
smoothly approaches the origin. 
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Figure 6: The radial probability density for the lS'1/2, 25*1/2: 35*1/2 f^nd 45*1/2 
states for a coupling of a = 0.35. The nodal pattern seen in figure|5lis beginning 
to get washed out as the wavefunction compresses around the horizon. 
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Figure 7: The radial probability density for the lS'1/2, 25*1/2: 351/2 ^nd 45i/2 
states for a coupling of a = 0.5. The pattern of nodes and dips seen in figures[Sl 
and|Blhas almost completely vanished, leaving a series of density profiles that 
lack structure. 
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Figure 8: The expectation value of r in units of GM/c^ for the S, P and D 
states. The broken horizontal line at x=2 shows the position of the horizon. 



in the obvious manner as 

/„°°dr(KP + |z.2P)- 

These are calculated via a straightforward Simpson's rule, and the results for 
the S, P and D orbitals are shown in figure|Hl We see that (r) decreases as the 
coupling increases. In the low-alpha regime, the expectation value follows the 
radius of the classical circular orbit, so (r) ex . As the coupling increases, 
and stable orbits become classically impossible, we find that (r) approaches, and 
moves within, the horizon. For higher a the bulk of the probability density lies 
inside the horizon, representing a short-lived state of a tightly-bound particle. 

While the low angular momentum orbitals are concentrated near the horizon, 
the orbitals with larger angular momentum still lie an appreciable distance out. 
As such, they adopt a form closer to the familiar hydrogen atom orbitals. A 
series of such orbitals are shown in figure El which shows the first excited mode 
for K values of 1, 2, 3 and 4. The coupling is again set to 0.5. As expected, the 
probability density is concentrated successively further from the hole. By the 
time we reach k = 3 (a classical radius of a; = 33) the wavefunction returns to 
the familiar hydrogen-like form. 



6 Decay rates 

So far we have concentrated on the real part of the energy, and the associated 
orbitals. But the fact that the black hole effective Hamiltonian is not Hermitian 
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Figure 9: The radial probability density for a range of angular momentum values 
with a coupling of a = 0.5. The first excited states are shown for k = 1,2, 3, 4. 
As K increases the orbitals are concentrated further from the source, and begin 
to resemble hydrogen atom wavefunctions. 
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Figure 10: The imaginary and real energies for the 15*1/2 state. The left hand 
plot shows (minus) the imaginary component of tlie energy as a function of the 
coupling strength. As expected, this increases as the orbits become more tightly 
bound. For comparison, the more complicated behaviour of the real part of the 
energy is shown on the right-hand side. 

implies that the energy is not real and the states have a finite half-life. As 
such the solutions could be viewed as representing resonance states as opposed 
to bound states. But for suitably large angular momenta the half lives can be 
pushed up as high as desired and the states will be extremely long lived. Such 
states are appropriate for a quantum description of a particle in a classically 
stable orbit some distance from the horizon. 

As argued above, the imaginary part of the energy will be negative, corre- 
sponding to a decay. The behaviour of this decay can be visualised in a number 
of ways. With E = lo ~ iv, the relevant quantity to study is 



In figure [Till we plot a as a function of coupling for the lS'1/2 state. For com- 
parison the real part of the energy is also plotted. The real energy falls to a 
minimum and starts increasing again as the orbits become unfavourably close, 
whereas the imaginary term simply increases monotonically. This as one would 
expect, as figure |S| showed that the orbits become increasingly tightly bound 
as a increases. As the coupling strength reaches 1, the imaginary component 
of the energy is of the order of 0.3 times the rest energy of the particle. This 
implies that the orbit should decay on the time-scale defined by the Compton 
frequency. These states are therefore extremely short lived, with a resonance 
width comparable to the orbital energy. 

In figure 1111 a is plotted for states with a range of angular momenta, k = 
1 • • • 5. The set of lowest-energy states (lS'1/2, 2P3/2,. . . ) is compared to the 
set of first-excited states (2S'i/2, 3P3/2,...). Both plots show the expected 
monotonic increase in a with coupling strength as the orbits become more tightly 
bound and a greater percentage of the wavefunction lies inside the horizon. The 
first-excited states are less tightly bound than the ground states, so have smaller 
decay rates. Below a threshold value of a, the imaginary energy is negligible. 
This threshold depends roughly linearly on k, and is the same for the lowest 
and first-excited states. Using the effective-potential model, we would expect 
decay to become dominant beyond the last value of a that allows stable circular 
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Figure 11: The imaginary energies of states with a range of angular momenta 
K = 1 • • • 5. The top plot shows the decay rates of the ground states, and the 
bottom plot shows the decay rates of the first-excited states, as functions of the 
coupling strength a. The positions of the minima in the real energy are marked 
with crosses. 
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orbits, a — n/\/V2 = 0.29k. The plot suggests that this model is reasonably 
valid. States with higher angular momentum can therefore be extremely stable, 
as the increase in n keeps the bulk of density away from the singularity. 

A classical argument can also be used to relate the high-a behaviour of the 
imaginary energy to the expectation value of wavefunction radius, by considering 
the proper time for radial infall. A massive particle starting at radius from 
rest would take proper time 



to reach the singularity. Conversely, the typical decay time for the wavefunction 
is ^ 

'^dccay — o V ' ' J 

If the decay time is similar to the infall time from the wavefunction expectation 
position (x), we would expect 

aacx(x)"^/^ (78) 

This model works well for the l>5'i/2 state, and we find aa cx {x) ^'^ in the 
high-a regime. The model requires some modifications for states with orbital 
angular momentum, as the infall time takes a more complicated form. 

With the decay rates now obtained, we can return to equation (|5Uf) to check 
the consistency of our method. For a number of states we computed the normal- 
ization integral and also extracted the behaviour of the state near the singularity. 
For all of these the imaginary component of the energy was consistent with equa- 
tion H5()|l . This confirms that the states are normalizable and represent genuine 
bound states. 



7 Discussion 

We have demonstrated the existence of a complicated spectrum of bound states 
for a quantum fermion in a black hole background. Each state represents a 
spatially-normalizablc solution to the Dirac equation in a Schwarzschild back- 
ground. The fact that time-separable solutions exist is simply established in 
one particular gauge, which casts the equation in a Hamiltonian-like form. A 
study of the behaviour of the wavefuntion under gauge transformations show 
that time-separability is a gauge-invariant feature. The spectrum itself is de- 
termined by boundary conditions applied at the horizon and at infinity. These 
alone are sufficient to imply the existence of an imaginary (decay) contribution 
to the energy. The physical explanation for this is provided by the singularity, 
which acts as a current sink. 

The qualitative features of the spectrum can be understood in terms of 
simple semi-classical models, but a full quantitative understanding only seems 
possible through a mixture of computational methods. The work in this paper 
can clearly be extended in a number of ways. We have only plotted the spectrum 
at low coupling strengths of a ~ 1, but astrophysical values can be far larger 
than this, with a ^ 10^^ for solar mass black holes. For larger a, the ground 
state will be one of high angular momentum. In this regime the spectrum will 
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be quite different to that of the hydrogen atom. One important question is 
precisely how great a binding energy can be achieved. In figure 0] we see that at 
around a = 5 we are achieving total energies of O.SSmc^, which is significantly 
lower the classical value of 0.94mc^. This suggests that more energy may be 
available in accretion processes than is traditionally thought. 

As well as increasing a, it would be of considerable interest to repeat this 
work for the case of a Kerr black hole. In this respect a useful start has been 
made in |15| . where the Kerr solution is written in a form which generalises 
the 'Newtonian' gauge employed in this paper. The calculations for the Kerr 
case are more complicated, however, because the angular separation constants 
are energy-dependent [3 IE]- There are also signs that the horizon structure 
of the Kerr solution will complicate the fairly straightforward picture presented 
here. The problem can be seen by analysing behaviour in a Reissner-Nordstrom 
background using the setup of this paper. For this case we find that the regular 
solutions at the outer horizon do not match onto regular solutions at the inner 
horizon. So quantum mechanics predicts that the probability density will pile up 
around the inner horizon in a similar manner to the classical picture. Behaviour 
of this type is inevitable, as the Reissner-Nordstrom singularity does not act as 
a sink, and the Hamiltonian is Hermitian. Since the current must still be inward 
pointing at the outer horizon, the probability density has to pile up somewhere. 
It seems likely that a similar picture holds for the Kerr solution, but detailed 
calculations are required to confirm this. 

The energy spectra presented in this paper raise a number of fundamental 
issues, which demonstrate the limitations in our current understanding of the 
interaction between gravity and quantum theory. It is unusual to obtain a decay 
law from quantum mechanics without some form of approximation. That we 
do so in the present case is a consequence of the fact that the system is open. 
States are allowed to decay onto the singularity, but no accompanying emission 
is considered. A complete treatment of the problem as a closed system would 
require a quantum theory of the singularity, and such a theory does not yet 
exist. 

The decay rates represent one feature of the quantum-mechanical description 
of the capture process. But, as well as decay, the quantum description of a 
particle falling onto the singularity of a black hole can involve a series of quantum 
jumps to lower energy orbits. This quantum description alters the physics of the 
process quite dramatically from the classical picture. As the particle undergoes 
a series of transitions we expect that it should radiate, which does not happen 
classically. Quite what form this radiation should take (electromagnetic, gravity 
waves?) is unclear. Also, as a transition takes place we should keep careful 
track of the evolution of the matter stress-energy tensor to tell us where the 
radiated energy is concentrated. A related problem this exposes is that we have 
not considered back reaction on the gravitational field, which could alter this 
picture. 

The quantum treatment of a particle in a gravitational field exhibits a curi- 
ous anti-parallelism with the electromagnetic case. In classical electrodynamics 
a charged particle in orbit around a point source should radiate, making atoms 
unstable. This problem is resolved by quantum mechanics, which predicts the 
existence of stable, non-radiating bound states. The reverse is true of gravi- 
tation. Classically, a particle can orbit a black hole in a geodesic outside the 
horizon, and such an orbit is stable. But quantum theory changes this, and 
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states that no totally stable orbits exist, due to the finite probability of the 
particle finding itself inside the horizon and ending on the singularity. While 
the time-scales involved in these decays may be of limited interest astrophys- 
ically, such processes are clearly of fundamental importance in understanding 
the interplay between quantum theory and gravitation. 

A final point to raise here is that the spectrum of real energies derived here 
has a mirror image of negative energy bound states. Each of these negative 
energy states also has a finite lifetime. If we model the vacuum in terms of a 
Dirac sea of filled negative energy states, we must include the bound states as 
well as the free continuum. It then follows that the vacuum itself is decaying 
— the black hole is sucking in the vacuum. Such a loss of negative energy 
states is seen as a creation of positive energy modes, which could contribute 
to Hawking radiation. This contribution appears to have been neglected in 
previous calculations, which concentrate only on the scattered states Q]. It 
is well known in calculations of the Lamb shift, for example, that ignoring the 
bound states in the calculation gives the wrong answer |E| ■ It would be of great 
interest to assess the contribution played by bound states to the gravitational 
analogue of this process. 
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